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FOREWORD 

This  research  paper  has  been  written  because  of  the  need  for 
explicit  characterization  of  smoothing  and  prediction  techniques 
for  many  applications.  While  the  methods  of  curve-fitting  by  least- 
squares  have  been  known  for  many  years,  it  is  frequently  difficult  to 
find  formulas  which  describe  the  errors  in  estimation  based  on  such 
methods . 

In  certain  applications,  the  central  problem  is  to  predict  the 
value  of  a  measured  quantity  which  exhibits  a  systematic  trend  which 
can  be  described  in  terms  of  a  low-order  polynomial.  In  other  appli¬ 
cations,  the  problem  is  to  provide  concurrent  estimates  of  the  true 
value  of  the  measured  quantity  and  its  rate  of  change.  In  both 
situations,  the  analytical  need  exists  for  measures  of  the  effects  of 
random  errors  and  their  interdependency  and  the  consequences  of 
systematic  errors  stemming  from  an  inadequate  model  of  the  underlying 
trend  This  paper  represents  an  initial  step  toward  fulfilling  this 
need.  It  was  motivated  by  problems  of  aiming  anti-aircraft  guns 
against  maneuvering  targets  and  achieving  precision  weapon  delivery 
by  tactical  aircraft.  It  is  hoped  that  the  results  presented  here 
will  find  other  applications  as  well. 


ABSTRACT 

Explicit  formulas  are  presented  for  estimating  position,  velocity, 
and  acceleration  in  low-order  polynomial  trends,  based  on  least -squares 
smoothing  of  sampled  data  accompanied  by  statistically  uncorrelated 
measurement  errors.  Formulas  are  also  given  for  interpolation  and 
prediction  of  position  and  velocity.  Expressions  for  the  variances 
and  covariances  of  consistent  position,  velocity,  and  acceleration 
estimates  are  given,  and  the  systematic  errors  accruing  from  use  of  a 
trend  estimation  basis  which  is  one  order  lower  than  the  actual  trend 
are  presented. 

One  interesting  result  is  that  the  normalized  correlation  between 
the  errors  in  an  estimate  of  current  position  and  those  in  an  estimate 
of  current  velocity  approaches  ,/%/ 2  when  the  number  of  measurements  in 
the  estimates  becomes  large. 

Finally,  the  problem  of  implementing  real-time  least -squares 
estimation  and  prediction  formulas  in  practical  systems  is  discussed. 

It  is  concluded  that  arithmetic  execution  time  requirements  can  be 
relaxed  by  generating  certain  sums  recursively,  and  that  data  storage 
requirements  can  frequently  be  eased  by  collapsing  the  raw  data  into 
short -term-average  samples. 
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I.  INTRODUCTION 

The  purpose  of  this  paper  is  to  document  some  ostensibly  well- 
known  results  of  elementary  statistics,  in  forms  which  will  hopefully 
prove  useful  to  persons  interested  in  the  smoothing  of  statistically 
stationary  data  sequences.  In  particular,  the  results  presented  de¬ 
scribe  the  properties  of  estimators  based  on  low-order  polynomial 
least-squares  fits. 

While  such  estimators  are  frequently  not  optimum  for  a  given 
setting,  they  are  realizable,  so  that  the  results  presented  here  con¬ 
stitute,  in  one  sense,  a  lower  bound  on  the  quality  of  estimation 
that  can  be  performed  in  a  given  situation.  Conversely,  if  a  least- 
squares  fit  does  happen  to  be  the  optimum  procedure  for  a  particular 
problem,  the  results  provide  upper  bounds  on  the  performance  of  more 
"economical"  procedures. 

In  what  follows,  it  is  assumed  that  the  data  originate  as  sampled 
values  of  a  continuous  well-defined  process,  x(t);  this  process  is 
assumed  to  be  representable  with  sufficient  precision  as  a  Taylor's 
series  in  the  independent  variable  t: 

x(t)  =  xQ  +  vQ(t-t0)  +  (l/2)a0(t-t0)2  +  (l/6)b0(t-t0)3  +  ...  (1) 

The  symbol  x  will  be  taken  to  represent  a  position  coordinate 
and  t  will  be  taken  to  represent  time,  so  that 

v(t)  =  dx/dt  (2) 


a(t)  =  d2x/dt2 


represent  the  velocity  and  acceleration  associated  with  the  coordinate 
x(t).  In  some  problems,  it  will  be  desirable  to  estimate  v  and  a  as 
well  as  x. 

The  approach  taken  in  this  paper  is  to  deal  with  the  underlying 
process  in  its  original  form,  as  given  by  Eq.  1,  rather  than  to  develop 
the  representation  of  Eq.  1  in  the  form  of  a  series  of  orthogonal  poly¬ 
nomials.  The  latter  approach  is  most  appropriate  when  the  desired  re¬ 
sult  is  simply  the  best  fitting  curve  for  estimation  and  prediction  of 
x(t),  and  it  is  used  extensively  in  statistical  literature,  e.g. ,  Ref. 
1,  pp.  186-191.  The  convenience  of  the  orthogonal-polynomial  approach 
stems  from  the  fact  that  the  problem  of  solving  the  simultaneous  equa¬ 
tions  for  the  coefficients  used  in  the  representation  is  trivial,  by 
virtue  of  the  orthogonality  property.  In  many  applications,  however, 
it  is  important  to  exhibit  the  consequences  of  the  terms  in  Eq.  1  as 
they  stand.  That  is,  position,  velocity,  and  acceleration  will  fre¬ 
quently  be  of  considerably  greater  significance  in  a  specific  situation 
than  the  coefficients  of  the  first,  second,  and  third  terms  in  the 
orthogonal -polynomial  representation.  It  is  for  this  reason  that  the 
more  cumbersome  direct  approach  has  been  adopted.  In  particular,  the 
direct  approach  facilitates  the  computation  of  variances  and  covari¬ 
ances  of  position,  velocity,  and  acceleration  estimates  and  the  deter¬ 
mination  of  the  consequences  of  systematic  errors  in  the  estimation 
model.  While  it  is  true  that  such  results  could  have  been  obtained 
from  the  orthogonal-polynomial  representation,  the  current  approach 
is  more  straightforward;  in  particular,  some  of  the  covariance  results 
are  more  easily  generalized  than  if  the  orthogonal-polynomial  approach 
had  been  adopted. 

The  true  sampled  values  are  unperturbed  values  of  x(t)  at  t  =  t 
+  kT,  where  k  is  an  integer  running  from  1  to  N  (N  is  therefore  the 
sample  size)  and  T  is  the  interval  between  adjacent  samples.  The  total 
period  over  which  the  data  are  observed  is 


r?*§g 

-y3£ 


Denoting  the  true  sampled  values  by  x^,  Eq.  1  gives 

\  =  x0  +  vokT  +  (l/2)aQ(kT)2  +  (l/6)b0(kT)3  +  .. 

For  simplicity  in  writing,  will  be  written  in  the  form 

2  3 

x.  =  x„  +  kr  +  k  q„  +  k  pA  +  ... 
k  o  o  no  ro 


with  the  correspondence 


(1/2)  aoT‘ 
(1/6)  bQT' 


and  so  forth. 

The  observed  data  will  be  denied  by  X^i  it  is  assumed  that  the 
errors  in  the  observed  data  manifest  themselves  as 


\  ~  \  +  ?k 


k  =  1,  ...  ,  N 


when  the  ?v’s  constitute  a  sequence  of  zero-mean,  uncorrelated,  iden- 
k  2 

tically  distributed  random  variables  with  variance  cx  • 

The  problem  at  hand  is  to  develop  means  for  estimating  the  true 
values  of  the  position,  or  velocity,  or  acceleration  (or  linear  func¬ 
tions  of  these  variables)  at  some  time  t,  which  will  be  expressed  in 
the  form 


t  =  t  +  XT 


.  -s 


(note  that  K  need  not  be  an  integer),  which  may  occur  prior  to,  during, 
or  subsequent  to,  the  period  of  observation.  Such  estimates  will  be 
denoted  by  the  symbols 


•  •  • 


for  the  estimated  position,  velocity,  acceleration,  etc.  at  the  time 
corresponding  to  the  choice  of  K. 

The  method  of  obtaining  formulas  for  these  estimates  will  be 
discussed  in  greater  detail  in  Sec.  II,  which  begins  with  an  assump¬ 
tion  about  the  nature  of  x(t).  This  assumption  deals  with  the  point 
at  which  the  Taylor’s  series  representation  of  Eq.  1  can  be  truncated 
without  sacrificing  the  usefulness  of  the  estimation  formula.  If  it 
is  assumed  that  the  values  of  x(t)  are  substantially  constant  over 
the  time  interval  for  which  the  estimation  formula  is  to  be  used, 
then  it  is  appropriate  to  employ  just  the  first  term  of  Eq.  1;  the 
estimation  formula  then  takes  the  form 

x^  =  xq  (zero  trend)  .  (9) 

If  it  is  assumed  that  x(t)  is  satisfactorily  represented  over  the 
estimation  interval  by  the  first  two  terms  of  Eq.  1,  then  the  estima- 
mation  formula  reads 

>L  =  x  +  r  K  (linear  trend)  .  (10) 

!\  O  U 

If  it  is  assumed  that  x(t)  is  satisfactorily  represented  over  the  es¬ 
timation  interval  by  the  first  three  terms  of  Eq.  1,  then  the  estima¬ 
tion  formula  reads 

x^  =  xQ  +  rQK  +  qQK2  (quadratic  trend)  ,  (11) 

and  so  forth.  For  reasons  that  will  become  apparent,  the  three  forms 
of  Eqs.  9,  10,  and  11  are  the  only  ones  that  will  be  considered  in 
this  paper. 
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The  least -squares  fitting  procedure  involves  choosing  the  param¬ 
eters  of  the  selected  estimation  formula  so  that  the  mean-square  dif- 

's  and  the  corresponding  x^’s  is  mini¬ 
mized.  Once  these  parameters  have  been  computed,  they  can  be  used  to 
provide  estimates  of  velocity,  acceleration,  and  so  forth.  Specifi¬ 
cally,  it  follows  that 

\  '  °>  *  0  (12) 
for  the  zero-trend  case; 

=  VT»  *K  =  0  (13) 

for  the  linear-trend  case;  and 

vK  =  (l/T)(ro+2qQK),  =  2qQ/T2  (14) 

for  the  quadratic -trend  case. 

The  remainder  of  the  paper  will  be  concerned  with  a  presentation 
of  formulas  for  the  estimation  parameters  (Sec.  II\;  an  assessment  of 
certain  statistical  characteristics  of  the  random  errors  in  estimation 
(Sec.  Ill);  a  discussion  of  the  systematic  errors  which  result  from 
an  inadequate  representation  in  the  smoothing  model;  and  finally,  a 
review  of  some  elementary  points  regarding  implementation  of  the  esti¬ 
mation  schemes  (Sec.  V). 

Some  clarification  of  the  notation  used  to  index  the  time  variable 
may  assist  in  interpreting  the  results  of  the  following  sections.  Three 
indices  are  employed: 

(1)  k,  which  is  an  integer,  is  used  to  index  the  times  at  which 
the  observed  data  were  taken; 

(2)  K,  which  is  a  continuous  variable,  is  used  to  index  the 
time  for  which  an  estimate  is  to  be  computed  and  was  chosen 
so  that  K  =  k  when  the  estimation  time  coincides  with  the 
time  at  which  the  kth  observed  datum  was  obtained; 


ference  between  the  observed 


(3)  M,  which  is  also  a  continuous  variable,  is  used  to  index  the 
time  for  which  an  estimate  is  to  be  computed,  and  was  chosen 
so  that  M  =  0  corresponds  to  the  midpoint  of  the  observation 
interval.  Because  the  observation  interval  extends  from 
K  =  1  to  K  =  N, 

M  =  K  -  (1/2)(N+1)  . 

All  three  indices  are  related  to  real  time  via  the  intersample  period 
T;  thus, 

t  =  t  +KT 
o 

=  tQ+MT  +  (1/2)(N+1)T  , 

where  t  is  the  value  of  t  at  one  intersample  period  prior  to  the 
time  of  observation  of  the  first  datum. 

In  some  applications  (prediction),  it  is  convenient  to  assume 
that  the  time  origin  coincides  with  the  time  at  which  the  last  (N^) 
datum  was  obtained.  Denoting  this  time  by  t', 


t'  =  t-t0-NT  , 


whence 


t '  -  (K-N)T 


t '  =  MT-T0/2  . 


i 


II.  ESTIMATION  FORMULAS 


A.  METHOD  OF  DERIVATION 

It  is  assumed  that  the  underlying  process  x(t)  can  be  represented 
by  a  truncated  Taylor’s  series,  e.g. ,  a  quadratic.  The  corresponding 
estimation  formula  is  similarly  truncated,  as  in  Eqs.  9,  10,  and  11. 
The  mean-squared  error  is  a  uniformly  weighted  average  of  the  squared 
difference  between  the  observed  data  and  the  value  yielded  by  the 
estimation  formula  at  the  corresponding  point  in  time: 


E  =  (\-V 


The  expression  for  E  that  results  when  the  explicit  form  of  the  esti¬ 
mator  xK  is  substituted  into  Eq.  15  is  then  separately  differentiated 
with  respect  to  each  of  the  (as  yet)  unknown  parameters  in  x^.  Because 
E  is  a  positive  quadratic  function  of  the  unknown  parameters,  equating 
the  derivatives  to  zero  yields  a  set  of  simultaneous  equations  which 
must  be  fulfilled  for  *  to  take  on  its  minimum  value.  These  equations 
are  then  solved  for  the  unknown  parameters.  In  solving  these  equa¬ 
tions,  the  following  identities*  are  helpful: 


]T)  k  =  (N/2)(N+1) 


(16) 


je ,  for  examf 


2,  pp.  7-8. 
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N 

k2  =  (N/6)(N+1)(  2N+1)  (17) 

k=l 

H 

]£)  k3  =  (N2/4)(N+1)2  (18) 

k=l 

N 

2  k4  =  (N/30)(N+1)(2N+1)(3N2+3N-1)  (19) 

k=l 


B.  RESULTS 
1.  Zero  Trend 


The  single  parameter  is  xq: 


ox> 

*-• 

• 

(20) 

From  Eq.  9,  this  is 
assumption , 

also  the  expression  for  x^.  Finally,  by 

o 

n 

«>* 

(21) 

o 

ii 

(22) 

2.  Linear  Trend 

The  two  parameters  are  xQ  and  rQ: 

N 

XQ  =  [  2/N(  N-l)]  2  [< 2N+1)  -  3k]  Xr 

J  k=l 

(23) 

8 


■-•<wa-i 


A 

r 


N 

[6/N(N2-l)j  ]|T  [2k  -  (N+l)]  ^ 
k=l 


(24) 


Substituting  these  results  into  the  position  estimation  formula, 

Eq.  10,  yields 

N 

^  =  [2/N(N2-l)]  2  [(N+1)(2N+1-3K)  +  3k(2K-N-l)]  (25) 

k=l 

For  the  special  case  K  =  N  (estimate  of  current  position),  Eq.  25 
simplifies  to 

N 

=  [2/N(N+l)j  2  [3k  -  (N+l)]  ^  .  (26) 

k=l 

The  velocity  estimate  is  simply 

iK  =  tyr  ■  (2?) 

Using  Eq.  4,  and  the  above  expression  for  rQ, 

N 

vK  =  [6/N(N+l)T0]  [2k  -  (N+l)]  \  ,  (28) 

k=l 


where  TQ  is  the  total  period  of  observation. 
-  By  assumption. 


(29) 


9 


<9  Jigs  «► 


3.  Quadratic  Trend 

The  three  parameters  are 


N 

{0  ■  Z  K  -  « »*l>k  +  < a)2+3N+2)]  \  <  30> 

k=l 


N 


•S  l 


A 

r  = 


°  N(N  -1)(N  -4)  " 


5 - S  [-30(N+l)k2  +  2( 2N+1)( 8N+ll)k 


-  3(  N+l)  (  N+2)  (  2N+1 ) J  \  (31) 


and 


N 

<3  - x-—-* - V  T6k2  -  6(N+l)k  +  (N+1)(N+2)1  X.  (32) 

0  N(N2-l)(N2-4)  “  L  J 

Substituting  these  results  into  the  position  estimation  formula, 

Eq.  11,  yields 


a  3 

*K  = 


N 

Y:  lo[(N+l)(N+2)  -  6X(N+1)  +  6K2]  k2 


N(N2-1)(N2-4) 

-2  ^3(N+l)(N+2)( 2N+1)  -  2K(2N+1)(8N+11)  +  30K2(N+1)]  k 

+( N+l) (N+2)  £(3N2+3N+2)  -  6K(2N+1)  +  10K2]  J  XR.  (33) 

For  the  special  case  X  =  N  (estimate  of  current  position),  Eq.  33 
simplifies  to 


N 

^N  =  Wm  WTTt  Y!  [10k2  -  2( 4N+3)k  +  (N+1)(N+2)J  XR  .  (34) 

k=l 
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Substituting  Eqs.  31  and  32  into  the  velocity  estimation  formula  of 
Eq.  14  and  invoking  Eq.  4,  there  is  obtained 


vv  - - ^ -  T1  3of2K  -  (N+l)l  k5 

*  N(N+1)(N  -4)TQ  "  L  J 


-2^30K(N+1)  -  (2N+1)(8N+11)J  k 
(N+l)(N+2)  [lOK  -  3(2N+1)J  | 


For  the  special  case  K  =  N  (estimate  of  current  velocity),  Eq.  35 
simplifies  to 


*  - - L- - V  f30(N-l)k2  -  2(14N2-11)  k 

N  N(N+l)(r-4)T0  “  L 


(N+1)(N+2)(4N-3)J 


The  acceleration  estimate  is  simply 


*K  *  V** 


Using  Eq.  4  in  Eq.  32, 


L  - - £0£N-12 -  V  [ek2  -  6(N+l)k  +  (N+1)(N+2)1  X.  .(37) 

*  N(N+1)(N  -4)T^  L  J  K 


4.  Interpretation  of  Results 


Each  of  the  estimation  formulas  involves  a  weighted  linear  aver¬ 
age  of  the  X^’s.  For  the  special  case  K  »  N  and  N»i,  the  relative 
weights  assigned  to  the  observed  data  for  current  position  estimation 
are  approximately: 


1 


zero  trend 


2  3(k/N)  -  lj  linear  trend 

3  10(k/N)2  -  8(k/N)  +  lj  quadratic  trend 

These  functions  are  depicted  graphically  in  Fig.  1.  It  can  be  seen 
that  as  the  order  of  the  estimation  formula  increases,  mere  weight 
is  placed  on  the  most  recent  data  (i.e. ,  those  points  for  k  approach 
ing  N). 

For  current  velocity  estimation  (K  =  N),  the  relative  weights 
for  N»!  are 

6  2(k/N)  -  lj  linear  crend 

12  15(k/N)2  -  14(k/N)  -s-  2]  quadratic  trend 

These  functions  are  depicted  graphically  in  Fig.  2.  It  is  seen  that 
in  the  quadratic  case,  the  more  recent  data  are  weighted  somewhat 
more  heavily  than  the  oldest  data,  and  that  the  extremes  are  weighted 
more  heavily  than  the  intermediate  points. 

Finally,  the  asymptotic  weighting  for  the  acceleration  estimate 
for  the  quadratic  case  is  shown  in  Fig.  3.  The  data  are  weighted 
symmetrically,  and  the  values  are  substantially  greater  than  for  the 
current  position  and  current  velocity  estimates. 
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FIGURE  1  Asymptotic  Weighting  for  Currant  Position  Estimation 
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III.  RANDOM  ERRORS 


A.  GENERAL  REMARXS 


The  estimation  formulas  presented  in  the  previous  section  are  all 
of  the  form 


A 

y 


N 


k=l 


(38) 


It  can  be  verified  that  these  formulas  have  the  following  property: 

if  the  model  (zero  trend,  linear,  or  quadratic)  for  the  underlying 

process  is  correct,  and  if  there  are  no  errors  in  the  observed  data, 

then  the  estimate  itself  is  correct.  In  this  section,  the  effect  of 

relaxing  the  second  condition  in  the  preceding  sentence  will  be  assessed. 

That  is,  it  will  be  presumed,  as  was  mentioned  in  the  Introduction, 

that  the  X^'s  are  in  error.  These  errors  are  assumed  to  be  additive, 

statistically  uncorrelated,  identically  distributed,  and  to  have  zero 

o 

means  and  (finite)  variance  o  .  Using  the  notation  of  Eq.  8,  the 

A  * 

estimate  y  is  given  by 


N 


N 


■  E  wk\ + L  % 


k=l 


k=l 


The  ensemble  average  of  y  is 


EP]  =Lvk  [h] 

k=l  k=l 
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Because  the  5^’s  are  zero-mean  random  variables,  the  second  sum  van* 
ishes.  In  view  of  the  preceding  remarks,  the  first  sum  is  just  the 
true  value  of  y,  so  that 

E  p]  *  y  • 

which  implies  in  turn  that  the  errors  in  y  have  zero  mean. 

The  next  step  is  to  compute  the  variance  of  the  errors  in  y.  It 
follows  from  the  foregoing  results  that  the  variance 


2(y)  =  efcy-y)2]  =  Ejj|>kS^j 

En  n 

E  E  w 

k=l  n=l 


N  N 


*  E  E  W  [?V?rJ 


k=l  n=l 


Because  the  errors  are  assumed  to  be  zero  mean  and  uncorrelated, 


M  = 


0  k  i  n 


ox  k  =  n 


so  that  Eq.  39  becomes 


The  task  of  computing  the  variance  of  the  estimation  formulas  is  thus 
reduced  to  one  of  identifying  the  weights  in  the  formula,  and  com¬ 
puting  the  sum  indicated  in  Eq.  40. 


~-******SW.»%5«*T-*  }■ 
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It  is  also  of  interest  to  compute  the  covariance  of  different 
errors.  Thus,  if  two  estimates  y,  given  by  Eq.  38,  and  z,  given  by 


.  n 

5  =  £  “k’Sc 


are  considered,  the  covariance  is  given  by 


Cov  f*y ,  zl  =  E  |Yy-y)  /z-zY| 


Following  the  same  procedure  as  for  the  variance  calculation, 


[*•  5]  =  °x2  E  WA 
k=l 


The  normalized  correlation  coefficient  of  the  errors  in  y  and  those 
in  z  is  given  by 


f')-S 


o(y)  a(£) 


B.  RESULTS 
1.  Zero  Trend 

The  variance  of  the  position  estimate  is 


'(\j  *  °2  ({o)  ■  °x7N 


By  assumption,  the  variances  of  the  velocity  and  acceleration  estimates 


are  zero. 
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2.  Linear  Trend 

The  variances  of  the  estimation  parameters  and  rQ  are 


The  covariance  of  $  and  r  is 

o  o 


Cov 


2°%  1 2N+1 ^ 

TJ  VTTT/ 

(43) 

12"x2 

N(N2-1) 

(44) 

,  2 
“6cx 

nTn-TT  * 

« 

The  variance  of  the  position  estimate  x^  can  be  computed  directly,  or 
by  noting  that  by  virtue  of  Eq.  10,  it  follows  that 

°2(*k)  ■  °2(*o)  ♦  *  Cov[4o-  fo]  ♦  • 


I 


N(N2-1) 


|(N+1)( 2N+1)  -  6K(N+1)  +  6K 


(45) 


For  the  special  case  K  =  N,  the  variance  of  the  estimate  of  current 
position  reduces  tc 

2 


°2(*n)  °  Hr  (tSt) 

The  variance  of  the  velocity  estimate,  v^,  is  simply 


(46) 
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-  - 

i  -  * . . 

Using  Eq.  4, 


W  ■£-(»)• 


The  normalized  correlation  coefficient  for  the  errors  in  the  current 
position  estimate  and  the  velocity  estimate  is 


p(v  *n)  = 


Quadratic  Trend 


h  k 


The  variances  of  the  estimation  parameters  xQ,  rQ,  and  are 


\  _  3°x  [~  3N24-3N4-2  1 

\°/“  N  L(N“1)(N-2)J 

2 

*/a\  _  12ax  f(2N+l)(8N+ll) 


(ir-l)(ir-4). 


The  covariances  are 


- 


fo.  *o]  ■ 

fo’%]  = 


180 0^ 


N(N2-l)(N2-4) 


180o  ( 2H+1) 
N(N-l)(N-2'r 


N(N-l)(N-2) 


-180ol 


N(N-1)(N  -4) 


The  variance  of  the  position  estimate  ^  is  given  by 


!(0  - 

\N  NfN-1 


N(N  -1)(N  -4) 


(N+1)(N+2)(3N  +3N+2) 


-  12( N+l)(N+2)( 2N+1)K  +  12( 7N2+15N+7)K2 


-  120(N+1)K3  +  6 OK  |  . 


Setting  K  =  N  in  Eq.  55,  the  variance  of  the  estimate  of  current 
position  reduces  to 


3o  2 

0  (*n)  =  N(N+l)(N+2)  ( 3N  ‘3N+2)  * 
The  variance  of  the  velocity  estimate  v^  is 

O  “ 


N2-4 


where  TQ  is  the  total  period  of  observation.  The  variance  of  the 
estimate  of  current  velocity  (K  =  N)  is 


-1X8N-1H 


?/*  \  _  12cx  tv-x\  hai-iK 

W  •  sj  v^/  [~ik 


The  normalized  correlation  coefficient  for  the  errors  in  the 
current  position  estimate  and  the  current  velocity  estimate  is 


'(v\)  ■ 


9(N-1)(N-2K  2N-1) 
( 8N-11)( 3N+3N-2) 


if  i 


>  : 


w- 


The  variance  of  the  acceleration  estimate 


* 

*X 


is 


(60) 


The  normalized  correlation  coefficient  for  the  errors  in  the  current 
position  estimate  and  the  acceleration  estimate  is 


P 


3( 3N2-3N+2) 


(61) 


The  normalized  correlation  coefficient  for  the  errors  in  the  current 
velocity  estimate  and  the  acceleration  estimate  is 


(62) 


C.  INTERPRETATION  OF  RESULTS:  LARGE-SAMPLE  LIMITS 

For  the  linear-trend  case,  the  variance  of  the  position  estimate 
x„  takes  on  a  minimum  value  for 


K 


o 


(63) 


For  this  value  of  K,  the  variance  is 

°X\) =  °*2/"  ’ 

which  is  the  same  as  the  variance  of  the  zero-trend  estimate  of  posi¬ 
tion.  This  result  implies  the  following: 


(1)  If  the  objective  of  the  estimation  procedure  is  to  estimate 
the  true  value  of  the  position  at  the  midpoint  of  the 


■*  *>  '■* 


'»Cw' 


observation  interval,  then  the  variances  yielded  by  the  zero- 
trend  formula  and  the  linear-trend  formulas  will  be  the  same.* 
(2)  if  the  zero-trend  hypothesis  is  true,  then  use  of  the  linear- 
trend  formula  for  position  estimation  will  yield  a  variance 
that  is  everywhere  greater  than  the  variance  that  would  be 
obtained  with  the  zero-trend  formula  with  the  exception  of 
the  single  point  given  by  Eq.  63.  This  point  is  in  conso¬ 
nance  with  the  general  statistical  principle  that  simply  in¬ 
creasing  the  complexity  of  the  hypothesized  trend  beyond  the 
necessary  level  will  only  degrade  the  quality  of  the  resulting 
estimate. 

With  regard  to  the  quadratic -trend  position  estimate,  the  variance 
takes  on  three  extreme  values  at  the  points 


_  N+l 
"T 


(64) 


and 


(65) 


The  variance  at  the  point  given  by  Eq.  64  is  a  local  maximum,  and 
has  the  value 


(66) 


The  variances  at  the  points  given  by  Eq.  65  are  minima  and  have  the 
value 


(67) 


Furthermore,  it  will  be  seen  in  the  next  section  that  the 
systematic  error  in  the  zero-trend  formula  vanishes  for  this 
case  when  a  strictly  linear  trend  is  present. 


r^>  ?  „  r 


The  variance  of  the  quadratic -trend  position  estimate  is  equal  to  the 
variance  of  the  linear-trend  position  estimate  at  the  points 


at  which  the  variances  are  simply 


*2(\)  -  VA 


For  the  quadratic-trend  case,  the  variance  of  the  velocity- 
estimation  error  takes  on  its  minimum  value  for 


K 

*o  2 


at  which  point  the  variance  is 


W  •=*-(»)• 


which  is  just  the  variance  of  the  velocity-estimate  error  for  the 
linear-trend  case. 

For  large  values  of  N  (i.e.  ,  NelOO),  the  expressions  for  the 
variances  and  correlation  coefficients  can  be  replaced  by  considerably 
simpler  approximations.  Of  particular  interest  are  the  nonralized 
correlation  coefficients,  in  the  linear-trend  case,  p(x^ ,  v^J 
approaches ^ 3/2  =  0.8660.  In  the  quadratic -trend  case,  the  following 
limiting  approximations  hold 


p(v  vN)  *  VV2  =  0.8660 

p^V  a  yfs/Z  =  0.7454 

P^N’  **)  *  =  0.9682 


The  fact  that  the  values  for  these  coefficients  are  rather  large  indi¬ 
cates  that  cross  terms  in  error  calculations  involving  (for  example) 
concurrent  position  and  velocity  errors  cannot  be  ignored. 

Figure  4  presents  graphs  of  the  large -sample  limit  of  the  standard 
deviation  of  the  position-estimation  error  normalized  so  that  the 
standard  deviation  for  the  zero-trend  case  is  unity. 

Figure  5  presents  graphs  of  the  large-sample  limit  of  the  standard 
deviation  of  the  velocity-estimation  error,  normalized  so  that  the 
standard  deviation  for  the  linear-trend  case  is  unity. 

These  graphs  indicate  that  the  errors  for  the  quadratic-trend 
case  grow  substantially  faster  than  for  the  linear-trend  case  when 
the  objective  is  to  predict  position  or  velocity  beyond  the  period  of 
observation  (K/N  >1).  In  practical  situations,  however,  it  is  fre¬ 
quently  possible  to  make  the  observation  period  substantially  longer 
in  the  quadratic-trend  case,  because  systematic  errors  in  the  predic¬ 
tion  equation  will  be  less  for  the  quadratic  case  than  they  would  be 
for  the  linear  case. 
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IV.  SYSTEMATIC  ERRORS 


A.  GENERAL  REMARKS 

In  many  practical  situations,  the  form  assumed  for  the  estimation 
formula  is  only  an  approximation  to  the  actual  variations  of  the  under¬ 
lying  process.  The  significance  of  higher-order  terms  in  the  Taylor’s 
series  representation  of  the  underlying  process  is  determined  by  the 
time  interval  over  which  the  estimation  rule  is  to  be  applied,  as  well 
as  the  magnitude  of  the  coefficients  of  these  terms.  For  the  purpose 
of  interpolation ,  i . e . ,  estimation  for  times  within  the  period  of  ob¬ 
servation,  the  significant  time  period  is  the  period  of  observation 
itself.  For  purposes  of  extrapolation  or  prediction,  the  prediction 
interval,  i.e. ,  the  elapsed  time  between  the  period  of  observation 
and  the  time  at  which  the  predicted  estimate  is  to  be  determined,  is 
usually  more  important. 

In  any  event,  the  presence  of  higher-order  terms  in  the  represen¬ 
tation  of  the  process  leads  to  errors  in  estimation  which  cannot  be 
removed  simply  by  improving  the  quality  of  the  observed  data  or  by 
taking  a  larger  sample  of  data.  To  distinguish  these  errors  from  the 
errors  which  accrue  from  imperfections  in  the  observed  data  (i.e. ,  the 
random  errors  discussed  in  the  previous  section),  they  are  referred 
to  as  systematic  errors. 

The  results  presented  here  are  intended  to  facilitate  a  prelim¬ 
inary  assessment  of  systematic  errors.  The  approach  taken  has  been  to 
assume  that  the  trend  of  the  true  data  contains  a  single  term  (of  the 
next -higher  order)  beyond  what  is  accounted  for  in  the  particular 
estimation  formula.  For  example,  a  cubic  trend  is  assumed  for  esti¬ 
mation  rules  based  on  the  quadratic  model.  As  such,  the  results  pro¬ 
vide  a  basis  for  establishing  necessary  (but  not  sufficient)  conditions 


i 


for  determining  the  validity  or  accuracy  of  a  particular  estimation 
procedure,  as  long  as  the  underlying  data  are  characterized  by  a  single, 
well-behaved  trend  which  is  representable  with  reasonable  accuracy  via 
Taylor’s  series.  It  should  be  noted  that  there  are  entirely  credible 
situations  in  which  the  underlying  process  is  not  characterizable  in 
this  manner;  for  example,  a  linear  trend  can  have  an  abrupt  change  in 
slope,  or  can  change  to  a  quadratic,  during  the  estimation  interval. 

The  results  presented  here  are  not  applicable  to  such  situations,  which 
necessitate  either  a  more  detailed  analytical  model,  or  computer  sim¬ 
ulation,  for  assessment  of  estimation  performance. 

B.  RESULTS 
1.  Zero  Trend 

The  position  estimate  is  a  constant  given  by  Eq.  20.  It  is 
assumed  that  the  underlying  process  exhibits  a  linear  trend,  so  that 
the  true  sample  values  are  of  the  form 


\  =  xo  +  kro  ' 


The  systematic  error  in  the  position  estimate  will  be  defined  as 


JL  =  Elx 


and  (for  this  case)  is  given  by 


A*K  =  <V0T/2)(N+1“2K)  » 


where  vQ  is  the  true  velocity. 

In  what  follows ,  it  is  of  interest  to  reference  the  estimation 
time  index  K  to  the  midpoint  of  the  period  of  observation.  This  is 
accomplished  by  setting 

M  a  K  -  .  (72 


The  quantity  KT  represents  the  time  for  which  the  estimate  is  to  be 
made,  relative  to  the  time  t„  mentioned  in  the  Introduction.  The 
quantity  MT,  in  turn,  represents  the  time  for  which  the  estimate  is  to 
be  made,  relative  to  the  midpoint  of  the  period  of  observation.  With 
this  shift ,  Eq.  71  becomes 


=  -MV0T  • 


2.  Linear  Trend 

The  position  estimate  is  given  by  Eq.  25;  it  is  assumed  that  the 
true  values  of  the  sampled  data  exhibit  a  quadratic  trend: 

2 

x.  =  x  +  kr  +  k  q  . 
k  o  o  o 

The  systematic  error  of  the  position  estimate  is  given  by 

=  -(a0T2/12)[6jc2  ‘  6K<N+1>  +  (N+1)(N+2)J  ,  (74) 

where  a  is  the  true  acceleration.  In  terms  of  the  shifted  index  M, 


aJK  “  -<a0T2/2)[M2  - 


The  velocity  estimate  is  given  by  Eq.  35.  The  systematic  error  is 


&vK  =  (a0T/2)(N+l  -  2K) 


Avv  =  Ma  T 
)v  o 


3.  Quadratic  Trend 


The  position  estimate  is  given  by  Eq.  33;  it  is  assumed  that  the 
true  values  of  the  sampled  data  exhibit  a  cubic  trend: 


\  ~  xo  +  kr0  +  k2qo  +  k\ 
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The  systematic  error  of  the  position  estimate  is  given  by 
AX*  =  (boT3/120)[(N+l)(N+2)(N+3)  -  2K( 6N2+15N+11) 

+  30K2(N+1)  -  20K3]  ,  (78) 

where  bQ  is  the  true  value  of  the  third  derivative  of  the  underlying 
process.  Expressed  in  terms  of  the  shifted  index  M, 

AXk  =  M(bQT3/120)( 3N2-7-20M2)  (79) 

The  velocity  estimate  is  given  by  Eq.  35;  the  systematic  error  is 


Avk  =  <boT2/60)[: 


30K2  -  30(N+1)K  +  (6N2+15N+ll)j 


or,  in  terms  of  the  shifted  index  M, 


Avk  =  (boT2/120)[60M2  -  (  3N2-7)]  . 


-a 
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V.  IMPLEMENTATION  CONSIDERATIONS 

In  some  situations  (particularly  those  involving  real-time  esti¬ 
mation),  direct  implementation  of  the  formulas  presented  in  Sec.  II 
can  lead  to  substantial  problems  from  the  standpoint  of  data  process** 
ing.  Specifically,  the  requirements  for  data  storage  and  computational 
speed  (or  parallel  arithmetic  hardware)  can  become  burdensome,  unless 
certain  simplifications  and  approximations  are  adopted. 

Consider,  for  example,  the  task  of  obtaining  smoothed  estimates 
of  the  current  range,  azimuth,  and  elevation  of  a  radar  target  from 
raw  data  provided  by  a  monopulse  tracking  radar.  Suppose  that  the 
radar  operates  at  a  pulse  repetition  frequency  of  10,000  pulses  per 
second,  and  that  smoothed  estimates  of  the  current  values  of  each  of 
the  three  coordinates  are  to  be  provided  at  a  rate  of  100  per  second, 
with  a  smoothing  interval  of  one  second  and  a  computational  lag  of  10 
milliseconds.  Under  these  conditions,  it  will  be  necessary  to  perform 
three  computations  of  the  form 


10,000 


A 

X  , 


*  T,  w 


k=l 


A 

one  hundred  times  per  second.  In  this  expression,  x^  is  the  j  esti¬ 
mate  of  the  value  of  the  coordinate  x  (i.e*  ,  range,  azimuth,  or  eleva¬ 
tion)  10  milliseconds  ago;  X,  .  is  the  k^  measurement  of  that  coor- 

K+3  a 

dinate  in  the  sample  to  be  used  in  computing  x. ;  and  W.  is  the  weight 
assigned  to  in  accordance  with  the  appropriate  estimation  formula 


from  Sec.  II.  For  example,  if  the  quadratic -trend  estimation  is  to  be 
employed,  Eq.  34  gives 

wk  ■  ~  ~  6»sjrxw*2) 

in  which  N  =  10,000. 

If  the  computations  were  to  be  carried  out  in  the  straightforward 
manner,  it  would  be  necessary  to  perform  30,000  multiplications  and 
additions  in  10  milliseconds,  which  would  require  the  completion  of 
a  single  multiply-and-add  cycle  in  one  microsecond,  assuming  parallel 
processing  of  the  three  sets  of  data.  Moreover,  such  computations 
would  entail  storage  and  shifting  of  30,000  words  of  raw  data.  While 
both  the  computational  speed  and  data  storage  requirements  are  within 
the  state  of  the  art ,  they  would  be  considered  burdensome  for  many 
applications. 

The  first  simplification  to  be  noted  stems  from  the  fact  that 
the  various  estimation  formulas  can  be  reduced  to  simple  linear  com¬ 
binations  of  the  following  sums: 


N 

so”> 

k=l 

N 

Sl(j)  =  2  kXk+j 
k=l 

N 

s2(0)  =Ek2xk+j 

k=l 

Updating  these  forms  can  be  accomplished  without  large-scale  computa¬ 
tions.  For  example,  if  estimates  are  to  be  provided  at  the  original 
data  rate,  the  following  simplifications  can  be  used: 
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S0(j+1)  =  s0(j)  +  -  Xj+1 


s^j^l)  =  s^j)  -  S0(j)  +  NXjj+j+1 

S2(j+1)  =  S2(j)  -  2Sx(j)  +  S0(j)  +  ^2X}i+j+1 

If  these  algorithms  were  to  be  employed  in  the  most  straightforward 
manner  for  the  computations  of  the  aforementioned  example,  it  would 
be  necessary  to  perform  approximately  3000  operations  in  the 
10-millisecond  computation  interval.  Assuming  parallel  operation  for 
the  three  sets  of  data,  100  microseconds  would  be  available  for  per¬ 
forming  the  ten  updating  operations  (seven  of  which  are  additions  or 
subtractions).  From  the  standpoint  of  computational  speed  requirements, 
this  rate  is  much  more  reasonable  than  the  previous  requirement,  but 
the  data  storage  requirement  is,  however,  slightly  greater. 

To  reduce  both  the  computational  burden  and  the  data  storage 
requirement,  it  may  be  permissible  to  replace  the  original  data  with 
short-term  averages.  Thus,  suppose  that  the  X's  are  averaged  over  a 
short  time  interval,  such  as  10  milliseconds,  and  that  the  resulting 
average  is  regarded  as  a  new  raw  datum  which  was  observed  at  the  mid¬ 
point  of  the  averaging  interval,  then 


t 


where  Nft  is  the  number  of  observations  used  in  computing  the  short-term 
average  (100  in  the  example  being  disucssed).  Equation  73  shows  that 
the  systematic  error  in  Xj  ,  due  to  a  linear  trend  in  the  raw  data, 
vanishes,  because  of  the  interpretation  of  X..*  as  having  been  taken 
at  the  midpoint  of  the  short-term  averaging  interval  (i.e.  ,  M  =  0). 

The  most  significant  systematic  error  will  therefore  arise  from 


31 


quadratic  and  higher  (even)  order  trends;  from  Eq.  75,  the  order  of 
magnitude  of  the  systematic  error  will  be 


=  aoV 


Thus,  if  the  systematic  error  is  to  be  less  than  one  foot,  the  accel- 

5  2 

eration  aQ  should  not  exceed  2.4x10  ft /sec  ,  or  approximately  7500  g. 
For  most  applications,  aQ  will  be  considerably  smaller  than  this 
figure,  so  that  the  systematic  errors  in  X^*  will  be  negligible. 

The  impact  of  this  on  the  data  processing  requirements  is  quite 
significant,  because  it  is  now  only  necessary  to  accumulate  three 
sums  of  the  form  of  X..*;  the  three  additions  per  100  microseconds  can 
be  done  serially  instead  of  in  parallel  (assuming  appropriate  buffer 
storage) . 
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111.  ftftONftOIIINft  MILITARY  ACTIVITY 


.••mier Expiicit  formulas  are  presented  for  estimating  position, 
velocity,  and  acceleration  in  low-order  polynomial  trends,  based  on 
least -squares  smoothing  of  sampled  data  accompanied  by  statistically 
uncorrelated  measurement  errors.  Formulas  are  also  given  for  in¬ 
terpolation  and  prediction  of  position  and  velocity.  Expressions 
for  the  variances  and  covariances  of  consistent  position,  velocity, 
and  acceleration  estimates  are  given,  and  the  systematic  errors 
accruing  from  use  of  a  trend  estimation  basis  which  is  one  order 
lower  than  the  actual  trend  are  presented. 

One  interesting  result  is  that  the  normalized  correlation  be¬ 
tween  the  errors  in  an  estimate  of  current  position  and  those  in 
an  estimate  of  current  velocity  approaches  Ji/2  when  the  number  of 
measurements  in  the  estimates  becomes  large. 

Finally,  the  problem  of  implementing  real-time  least-squares 
estimation  and  prediction  formulas  in  practical  systems  is  dis¬ 
cussed.  It  is  concluded  that  arithmetic  execution  time  requirements 
can  be  relaxed  by  generating  certain  sums  recursively,  and  that 
data  storage  requirements  can  frequently  be  eased  by  collapsing  the 
raw  data  into  short -term-average  samples. 
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